﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using alglib;

namespace Engine.Mathematics
{
    /// <summary>
    /// Class for calculating conductance of lattice for non-linear elements
    /// </summary>
    public class NonLinearGMRES
    {
        private const int m = 100;

        /**
        * Givens rotations for the QR factorization
        */
        private GivensRotation[] rotation;

        private DenseMatrix _A;

        private Vector _B;

        private Matrix _H;

        private Vector _R;

        private IList<Vector> _V;

        /// <summary>
        /// Low determining behaviour of elements
        /// </summary>
        private NonLinearConductance _conductance;

        /// <summary>
        /// Bonds for this lattice
        /// </summary>
        private IList<Bond> _bonds;

        /// <summary>
        /// Points for this lattice
        /// </summary>
        private IList<Point> _points;

        /**When norm of vector _R is lower that 
         * DELTA we stop iterating throw this algorithm. */
        private const double DELTA = 1E-10;

        /** Initial approximation. */
        private Vector _x0;

        public NonLinearGMRES(NonLinearConductance conductance, IList<Point> points, IList<Bond> bonds)
        {
            this._points = points;
            this._bonds = bonds;
            this._conductance = conductance;
            //initial value of x
            _x0 = new Vector(_points.Count);
            CalculateLatticeMatrix();


            _V = new List<Vector>(_B.Size);
        }


        /// <summary>
        /// Calculate lattice new matrix
        /// based on current in-progress solution
        /// </summary>
        private void CalculateLatticeMatrix()
        {
            //we use this to track about entities that interested us in
            //calculation conductance
            this._A = new DenseMatrix(_points.Count, _points.Count);
            this._B = new Vector(_points.Count);
            foreach (Point point in _points)
            {
                int mainIndex = _points.IndexOf(point);
                //then find all relation with other points
                IList<Bond> bonds = point.BondList;
                double coefficientForPoint = 0;
                double partB = 0;
                foreach (Bond bond in bonds)
                {
                    Point nextPoint = point.Navigate(bond);
                    if (nextPoint != null)
                    {
                        int nextPointIndex = _points.IndexOf(nextPoint);
                        double bondConductance = BondConductance(mainIndex, nextPointIndex, nextPoint);
                        coefficientForPoint += bondConductance;
                        if (nextPoint.Potential == null)
                        {
                            _A[mainIndex, nextPointIndex] = 0d - bondConductance;
                        }
                        else
                        {
                            partB += bondConductance * (double)nextPoint.Potential;
                        }
                    }
                }
                _A[mainIndex, mainIndex] = coefficientForPoint;
                _B[mainIndex] = partB;
            }
        }

        /// <summary>
        /// Calculate bond conductance
        /// </summary>
        /// <param name="bond">bond for which we calculate conductance</param>
        /// <param name="nextPoint">next point to which connected this bond</param>
        /// <returns>conductance of bond</returns>
        private double BondConductance(int mainIndex, int nextPointIndex, Point nextPoint)
        {
            double voltage = 0;
            if (nextPointIndex == -1)
            {
                voltage = (double)nextPoint.Potential - _x0[mainIndex];
            }
            else
            {
                voltage = _x0[nextPointIndex] - _x0[mainIndex];
            }
            if (voltage != 0)
            {
                double conductance = _conductance(voltage) / voltage;
                return conductance;
            }
            else
            {
                return 1;
            }
        }

        /// <summary>
        /// Finds conductance of lattice with non-linear elements.
        /// </summary>
        /// <returns>Conductance of lattice with non-linear elements</returns>
        public Vector Solve()
        {
            //then we start our iterative algorithm
            //iterate to next stage of our algoritm
            //initial conditions.
            rotation = new GivensRotation[m + 1];
            _R = _B - (_A * _x0);
            Vector v0 = _R / (_R.Norm);
            _V.Add(v0);
            while (_R.Norm2 > (_x0.Norm > DELTA ? _x0.Norm * DELTA : DELTA))
            {
                Vector s = new Vector(m + 1);
                s[0] = _R.Norm;
                //create new matrix for storing 
                _H = new Matrix(m + 1, m);
                for (int j = 0; j < m; j++)
                {
                    Vector vTemp = _A * _V[j];
                    for (int i = 0; i <= j; i++)
                    {
                        _H[i, j] = (_A * _V[j]) * _V[i];
                        vTemp = vTemp - (_V[i] * _H[i, j]);
                    }

                    _H[j + 1, j] = vTemp.Norm;
                    Vector vNext = vTemp / vTemp.Norm;

                    // QR factorization of H using Givens rotations
                    for (int k = 0; k < j; ++k)
                        rotation[k].apply(_H, j, k, k + 1);

                    rotation[j] = new GivensRotation(_H[j, j], _H[j + 1, j]);
                    rotation[j].apply(_H, j, j, j + 1);
                    rotation[j].apply(s, j, j + 1);
                    _V.Add(vNext);
                }

                int info = -1;
                double[] solution = new double[m];
                densesolver.densesolverreport report = new densesolver.densesolverreport();
                double[,] RDataN = _H.Data;
                double[] RightPartOfEquationData = s.Data;
                densesolver.rmatrixsolve(ref RDataN, m, ref RightPartOfEquationData,
                    ref info, ref report, ref solution);
                if (info != 1)
                {
                    throw new Exception("Exception occured while solving this equation.");
                }
                Vector y = new Vector(solution);
                //them we need to estimate x for the next cycle.
                for (int i = 0; i < m; i++)
                {
                    _x0 = _x0 + (_V[i] * y[i]);
                }
                //prepare for the next cycle
                //iterate to next stage of our algoritm
                //initial conditions.
                _V.Clear();

                //recalculate a matrix
                CalculateLatticeMatrix();

                _R = _B - (_A * _x0);
                Vector v = _R / (_R.Norm);
                _V.Add(v);
            }
            return _x0;
        }

        /**check were we canceled our operation */
        public void Cancel()
        {
            //todo this method will check if user is canceled this operation.
            throw new NotImplementedException();
        }

    }
}
